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Complex low-temperature ordered states in chiral magnets are typically governed by a competition 
between multiple magnetic interactions. The chiral-lattice multiferroic Cu20Se03 became the first 
insulating helimagnetic material in which a long-range order of topologically stable spin vortices 
known as skyrmions was established. Here we employ state-of-the-art inelastic neutron scattering 
(INS) to comprehend the full three-dimensional spin excitation spectrum of Cu20Se03 over a broad 
range of energies. Distinct types of high- and low-energy dispersive magnon modes separated by an 
extensive energy gap are observed in excellent agreement with the previously suggested microscopic 
theory based on a model of entangled CU4 tetrahedra. The comparison of our INS data with model 
spin-dynamical calculations based on these theoretical proposals enables an accurate quantitative 
verification of the fundamental magnetic interactions in Cu20Se03 that are essential for understand¬ 
ing its abundant low-temperature magnetically ordered phases. 


Introduction. Chiral magnets form a broad class of mag¬ 
netic materials in which long-range dipole interactions, mag¬ 
netic frustration, or relativistic Dzyaloshinskii-Moriya (DM) 
interactions twist the initially ferromagnetic spin arrange¬ 
ment, thus leading to the formation of noncollinear incom¬ 
mensurate helical structures with a broken space-inversion 
symmetry^“"^. Helical magnetic states occur in a very broad 
range of compounds including metals and alloys, semicon¬ 
ductors, and multiferroics. From the fundamental point of 
view, the signihcance of such materials is due to the richness 
of their possible magnetic structures as compared to ordinary 
(commensurate and cohinear) magnets. Perhaps the most 
prominent example is given by the topologically nontrivial 
long-range ordered states known as skyrmion lattices^”^. 

The interest in the multiferroic ferrimagnet Cu20Se03 
with a chiral crystal structure has been intensified in re¬ 
cent years after the discovery of a skyrmion-lattice phase in 
this system^, thus making it the first known insulator that 
exhibits a skyrmion order. It was shown soon afterwards 
that this skyrmion lattice can be manipulated by the ap¬ 
plication of either magnetic^ or electric^^ field, and more 


recently also by chemical doping^which indicates a deli¬ 
cate balance of the magnetic exchange and DM interactions 
with the spin anisotropy effects. Understanding the complex 
magnetic phase diagram of Cu20Se03 therefore requires 
detailed knowledge of the spin Hamiltonian with precise 
quantitative estimates of all interaction parameters, which 
can be obtained from measurements of the spin-excitation 
spectrum. Microscopic theoretical models that were re¬ 
cently proposed for the description of spin arrangements in 
Cu20Se03 include five magnetic exchange integrals and hve 
anisotropic DM couplings between neighbouring S = | cop¬ 
per spins^^”^^, yet their experimental verihcation was so far 
limited to thermodynamic data^"^, terahertz electron spin res¬ 
onance (ESR)^^, far-infrared^^ and Raman spectroscopy^^ 
that can only probe zone-center excitations in reciprocal 
space. Here we present the results of inelastic neutron scat¬ 
tering (INS) measurements that reveal the complete picture 
of magnetic excitations in Cu20Se03 accessible to modern 
neutron spectroscopy over the whole Brillouin zone and over 
a broad range of energies and demonstrate good quantitative 
agreement with spin-dynamical model calculations both in 



Figure 1 | Simplified structure of the magnetic 
Cu sublattice in Cu20Se03 and Brillouin zone 
unfolding, a, The sublattice of Cu atoms approx¬ 
imated by an fee arrangement of identical tetra¬ 
hedral CU4 clusters, the corresponding unit cell 
shown with grey lines. Introduction of imaginary 
tetrahedra (faded color) within the voids allows 
to introduce a twice smaller structural unit cell 
(dashed red lines) that corresponds to a half-filled 
fee lattice of Cu atoms, b, Brillouin zones that 
correspond to the original crystallographic simple- 
cubic unit cell with the lattice parameter a (inner 
grey cube), the fee unit cell with lattice parameter 
a (smaller truncated octahedron, blue lines), and 
the fee unit cell with lattice parameter a/2 (larger 
truncated octahedron, red lines). 
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Figure 2 | Magnon excitations in Cu20Se03 as measured by triple-axis neutron spectroscopy in the low-temperature spin-spiral 
state (T = 5 K). a, Dispersion of the low-energy ferromagnon branch along the (OOH) direction mapped out with the thermal-neutron 
spectrometer PUMA. Data points represent peak positions that resulted from fitting the energy profiles for every momentum, b, The lower 
part of the same dispersion within the region shown by a black rectangle in panel (a), measured with higher resolution at the cold-neutron 
spectrometer PANDA. The magnon corresponds to the higher-energy broadened peak (squares), the sharp weaker peaks below it originate 
from spurious Bragg scattering. The peak positions from the thermal-neutron data in panel (a) are also shown here for comparison (grey 
circles), c, Raw energy scans across the upper magnon bands measured at several high-symmetry points using the PUMA spectrometer. 
Here the spectra are grouped and marked according to the small (original crystallographic) cubic Brillouin zone, d, Individual energy cuts 
from the cold-neutron dataset in panel (b) for different momenta along (2 2 2 - 5 ). The fits shown with solid lines neglect the spurious 
peaks at the low-energy side of the spectrum. 


terms of magnon dispersions and dynamical structure fac¬ 
tors. 

Crystal structure and Brillouin zone unfolding. The cubic 
copper(II)-oxoselenite Cu20Se03 crystallizes in a complex 
chiral structure with 16 formula units per unit cell (space 
group P 2 i 3 , lattice constant a = 8.925 A) This structure 
has been visualized, for instance, in Refs. 8 or 14 . For the 
discussion of magnetic properties relevant for our study, it 
is useful to consider only the magnetic sublattice of copper 
ions, which can be approximated by a face-centered-cubic 
(fee) lattice of identical CU4 tetrahedra, as shown in Fig. la. 
This reduces the volume of the primitive unit cell fourfold 


with respect to the original simple-cubic structure, resulting 
in only four Cu atoms per primitive cell. In reciprocal space, 
this simplified lattice possesses a larger “unfolded” Brillouin 
zone with a volume of 4 • ( 27 i/a)^ as shown in Fig. lb with 
blue lines. 

Existing theoretical models^^"^^ suggested that individual 
CU4 tetrahedra represent essential magnetic building blocks 
of Cu20Se03. Due to the structural inequivalence of the 
copper sites within every tetrahedron, its magnetic ground 
state orients one of the four Cu^^ spins antiparallel to three 
others that are coupled ferromagnetically, resulting in a state 
with the total spin S = 1 . The interactions between neigh¬ 
boring tetrahedral clusters are at least 2.5 times weaker than 
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intratetrahedron couplings^^”^^ and lead to a ferromagnetic 
arrangement of their total spins below the Curie tempera¬ 
ture, Tq = 57 K. In addition, weak DM interactions twist 
the resulting ferrimagnetic state into an incommensurate 
helical spin structure propagating along the (001) direction 
with a pitch = 63 nm, that is 100 times larger than the 
distance between nearest-neighbor tetrahedra, which cor¬ 
responds to the propagation vector = 0.010 derived 
from small-angle neutron scattering (SANS) data^. 

As the next step, we note that the same structure can be 
represented as a half-filled fee lattice of individual Cu atoms. 
Indeed, if one supplements the lattice with an equal number 
of imaginary CU4 tetrahedra by placing them within the voids 
of the original structure, as shown in Fig. la, a twice smaller 
fee unit cell with lattice parameter a/2 can be introduced 
(red dashed lines). Its primitive cell contains on average 
one half of a copper atom and corresponds to the large “un¬ 
folded” Brillouin zone with a volume of 32 • ( 27 i/a)^ that 
is shown with red lines in Fig. lb. As will be seen in the 
following, the dynamical structure factor representing the 
distribution of magnon intensities in Cu20Se03 inherits the 
symmetry of this large Brillouin zone. The described un¬ 
folding procedure is therefore useful for reconstructing the 
irreducible reciprocal-space volume for the presentation of 
INS data. Throughout this paper, we will denote momentum- 
space coordinates in reciprocal lattice units corresponding to 
the crystallographic simple-cubic unit cell (1 r. l.u. = 2 n/a), 
whereas high-symmetry points will be marked in accordance 
to the large unfolded Brillouin zone as explained in Fig. lb. 
Equivalent points from higher Brillouin zones will be marked 
with a prime. 

Triple-axis neutron spectroscopy. We start the presenta¬ 
tion of our experimental results with the triple-axis spec¬ 
troscopy (TAS) data shown in Fig. 2 . At low energies, the 
spectrum is dominated by an intense Goldstone mode with 
a parabolic dispersion, which emanates from the r\222) 
wave vector and has the highest intensity at this point, as 
can be best seen from the thermal-neutron data in Fig. 2 a. 
Because this point corresponds to the center of the large un¬ 
folded Brillouin zone, we can associate this low-energy mode 
with a ferromagnetic spin-wave branch anticipated from the 
theory for the collinear magnetically ordered state^^. Upon 
moving away from the zone center along the (001) direction, 
its dispersion reaches a maximum of about 12 meV, while its 
intensity is reduced continuously towards the unfolded-zone 
boundary. Much weaker replicas of the same ferromagnon 
branch can be also recognized at other points with integer 
coordinates, like (221) or (220), as they coincide with zone 
centers of the original crystallographic Brillouin zone but not 
the unfolded one. We have fitted individual energy scans in 
Fig. 2 a with Gaussian peak profiles to extract the measured 
magnon dispersion shown as black dots. 

For a closer examination of the low-energy part of the 
ferromagnon branch around the r^(222) point, we also per¬ 
formed TAS measurements at a cold-neutron spectrometer 
providing higher energy resolution, as shown in Figs. 2 b and 
2 d. Here we see a parabolically dispersing ferromagnetic 
spin wave branch (squares) in good agreement with the dis¬ 
persion obtained from the thermal-neutron measurements 
(grey circles). The data are additionally contaminated with 
spurious Bragg scattering that appears as a sharp straight line 
below the magnon peak in Fig. 2 b or as a strong low-energy 


peak in 2 d. While fitting the experimental data to extract 
the magnon dispersion (solid lines in 2d), this spurious peak 
had to be masked out. As a result, the ferromagnon mode 
appears to be gapless within our experimental resolution, 
which agrees with the earlier microwave absorption mea¬ 
surements, where the spin gap value of only 3 GHz (12 jUeV) 
has been reported^^. By fitting the experimental dispersion 
at low energies with a parabola, 1 t.cl> = Dq^, where q is the 
momentum transfer along (001) measured relative to the 
( 222 ) structural Bragg reflection, we evaluated the spin- 
wave “stiffness” D = 52 . 6 meV-A^, which is comparable to 
that in the prototypical metallic skyrmion compound MnSi^\ 

The magnetic Goldstone mode is significantly broadened 
in both energy and momentum and has a large intrinsic 
width that is considerably exceeding the instrumental reso¬ 
lution of ~ 0 . 14 meV This can be naturally explained as a 
result of the incommensurability of the spin-spiral structure 
and its deviation from the collinear ferrimagnetic order. In 
the helimagnetic state, the low-energy magnetic excitations 
are expected to split into so-called helimagnons emanat¬ 
ing from the 1 ^^-, 2 ^^-, and higher-order magnetic Bragg 
reflections^^"^"^. The incommensurability parameter is 
too small in our case to be resolved with a conventional 
triple-axis spectrometer. The typical helimagnon energy 
scale in Cu20Se03 can be estimated as Dk^ ^ 10.5 [leV, that 
is an order of magnitude lower than in MnSi and also lies 
well below the energy resolution of the instrument. As a 
result, multiple unresolved helimagnon branches merge into 
a single broadened peak as seen in Fig. 2 d. 

According to the theoretical calculations^^, a higher- 
energy band of dispersive magnon excitations is expected 
between approximately 25 and 40 meV In Fig. 2 c, we show 
several spectra measured in this energy range at a number 
of high-symmetry points that reveal at least three distinct 
peaks. Here we have grouped the spectra according to their 
location in the original crystallographic Brillouin zone (grey 
cube in Fig. lb) to emphasize that the peak intensities are 
strongly affected by the dynamical structure factors, whereas 
their positions in energy in different Brillouin zones remain 
essentially the same for equivalent wave vectors. The theory 
anticipates that these high-energy modes represent an intri¬ 
cate tangle of multiple magnon branches, so that several of 
them may contribute to every experimentally observed peak. 
A direct comparison with the existing spin-dynamical calcu¬ 
lations therefore appears difficult with the limited TAS data 
and requires a complete mapping of the momentum-energy 
space using time-of-flight (TOF) neutron spectroscopy that 
we present below. 

Time-of-flight neutron spectroscopy. The benefit of TOF 
neutron scattering is that it provides access to the whole 
4 -dimensional energy-momentum space in a single 

measurement, which is particularly useful for mapping out 
complex magnetic dispersions that persist over the whole 
Brillouin zone. The data can be then analyzed to extract any 
arbitrary one- or two-dimensional cut from this dataset as 
we describe in the Methods. Fig. 3 presents several typical 
energy-momentum cuts taken along different high-symmetry 
lines parallel to the (001), (110), and (111) crystallographic 
axes. One clearly sees both the low-energy ferromagnon 
band that is most intense around r'{222) and the inter¬ 
twined higher-energy dispersive magnon modes between 
25 and 40 meV whose intensity anticorrelates with that of 
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Figure 3 | Representative energy-momentum cuts through 
the TOF data along the following high-symmetry directions: 

a, (IIL); b, (22L); c, (HHl); d, (HH2); e, (HHH); f, (12L). All 
the presented data were measured at the base temperature of 4.5 K. 
Black and white arrows mark two high-energy optical phonon 
modes (see text). 

the low-energy modes: It vanishes near r^(222) yet is maxi¬ 
mized near theX^(220) point at the unfolded-zone boundary, 
where the low-energy modes are weak. The low- and high- 
energy modes are separated by a broad energy gap in the 
range ~13-25meV In addition, at much higher energies 
one can see two relatively weak flat modes centered around 
49 and 61 meV (marked with arrows in Fig. 3a). Although 
the theory predicts a flat magnetic mode around 54 meV^^, 
the fact that their intensity monotonically increases towards 
higher |Q| speaks rather in favor of optical phonons. It is 
likely that the 54 meV magnetic mode is far too weak com¬ 
pared to the phonon peaks and therefore cannot be clearly 
seen in our data. As the data in Fig. 3 extend over several 
equivalent Brillouin zones in momentum space, one can see 
an increase of the nonmagnetic background and a simulta¬ 
neous decrease of the magnetic signal towards higher |Q|, 


resulting in a rapid reduction of the signal-to-background 
ratio. Therefore, in the following we restrict our analysis 
mainly to the irreducible part of the reciprocal space corre¬ 
sponding to shortest possible wave vectors between r(OOO) 
and r(222). 

It is interesting to note how the hierarchical structure 
of Brillouin zones introduced in Fig. 1 is directly reflected 
in the magnon intensities. Along the (HHH) direction in 
Fig. 3e, the strongest spin-wave modes are seen around cen¬ 
ters of the large unfolded zones, r(OOO) and r'(222), while 
a somewhat weaker mode appears at the L(lll) point that 
coincides with the zone center for the smaller fee Brillouin 
zone (shown in blue in Fig. lb). Other points with integer co¬ 
ordinates that are centers of the smallest simple-cubic zones 
but of none of the larger ones, such as (110) or (001), con¬ 
tain only much weaker replica bands that are barely seen in 
the data (e.g. Fig. 3a,c or Fig. 2a). This leads to the appear¬ 
ance of rather unusual features in certain cuts through the 
spectrum, such as the one at the (001) point in the middle 
of Fig. 3c, which visually resembles the bottom part of the 
famous hourglass-like dispersion in some transition-metal 
oxides^^”^^. However, here its intense top part stems from 
the nearest unfolded zone center r(OOO) that lies outside 
of the image plane, while the weaker downward-dispersing 
branches originate from the L points at (111) and (111), 
producing such a peculiar shape. The situation is even more 
complicated at higher energies, where the strongest well- 
resolved downward-pointing modes with a bottom around 
26 meV are located at W points of the large unfolded zone, 
such as (102), (210) or (320), with weaker replicas shifted 
by the (111) or equivalent vector. Interestingly, without 
Brillouin zone unfolding such a particularity of the W points 
would not at all be obvious without explicit calculations of 
the dynamical structure factors, as they are all equivalent 
to the zone center in the original crystallographic Brillouin 
zone notation. 

The pattern of INS intensity in momentum space is visual¬ 
ized by the constant-energy slices presented in Fig. 4. All of 
them are oriented parallel to the (HHL) plane, which was 
horizontal in our experimental geometry. The low-energy 
cut integrated around 8 ± 1 meV (Fig. 4a) passes through 
the origin and contains both the (222) and (111) wave vec¬ 
tors, where rings of scattering from the stronger and weaker 
magnetic Goldstone modes, respectively, can be clearly seen. 
The next cut at 28 ± 1 meV in Fig. 4b is shifted by the (||0) 
vector in such a way that it crosses the previously mentioned 
downward-pointing high-energy modes, which appear as 
small circles around W points. Finally, the last two panels 
(Fig. 4c,d) show integrated intensity in the (HHL) plane 
at energies corresponding to the middle (32.5 ±3.5 meV) 
and top (38 ± 2meV) of the upper magnon band, respec¬ 
tively. While individual magnon dispersions are intertwined 
at these energies and cannot be resolved separately, in both 
panels we observe deep minima of intensity at the r\222), 
r''(400), and equivalent wave vectors. They correspond 
to the suppression of the dynamical structure factor at the 
center of the unfolded Brillouin zone, where the low-energy 
ferromagnon mode, on the contrary, is most intense. 

To get a more complete picture of magnetic excitations 
in Cu 20 Se 03 and to compare them with the results of spin- 
dynamical calculations, we have symmetrized the TOF data 
(see Methods) and assembled energy-momentum profiles 
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Figure 4 | Constant-energy slices through the TOF data parallel to the (HHL) plane, a, Along the (HHL) plane at Hco = 8 ± 1 meV 
cutting through the low-energy magnon band with intense rings of scattering from the ferromagnon modes around F points; b, Along the 
offset (H—I H-l-1 L) plane at Hco = 28 ± 1 meV through the bottom part of the upper magnon band, showing rings of intensity around W 
points; c, d, Along the (HHL) plane through the upper magnon band at Hco = 32.5 ± 3.5 meV and 38 ± 2 meX respectively. In all panels, 
the data were integrated within ±0.15 r. l.u. from the indicated plane along the (110) direction that is orthogonal to the plane of the 
figure. Corresponding energy integration ranges given above each panel are also shown in Fig. 5a with vertical bars. 


along a polygonal path involving all high-symmetry direc¬ 
tions in Q-space into a single map presented in Fig. 5a. Due 
to the kinematic constraints, the data in the first Brillouin 
zone are always limited to lower energies as can be seen 
in Fig. 3, but have the best signal-to-noise ratio. Therefore, 
for the best result we also underlayed data from higher Bril¬ 
louin zones to show the higher-energy part of the spectrum 
at equivalent wave vectors in the unfolded notation. 
Spin-dynamical calculations. Based on previous theoret¬ 
ical and experimental results^^”^^, we used a tetrahedron- 
factorized multi-boson method to calculate the magnon spec¬ 
trum for the collinear ferrimagnetically ordered state. We 
consider the Heisenberg-type Hamiltonian: 

^ = J]jijSi-Sj, ( 1 ) 

ij 

where Jij denote five different kinds of exchange couplings 
between sites i and j that are explained in Fig. 5c. There 
are two strong couplings, a ferromagnetic (J™) and an an¬ 
tiferromagnetic (J^) one, residing on the CU 4 tetrahedra 
that form the elementary units of our tetrahedra-factorized 
multiboson theory. Furthermore, there are weak ferromag¬ 
netic (J™) and antiferromagnetic (J^^) first-neighbor inter¬ 
actions connecting the tetrahedra. The last coupling is an 
antiferromagnetic longer-range interaction (J0..0) connect¬ 
ing Cui and CU 2 across the diagonals of hexagonal loops 
formed by alternating Cu^ and CU 2 sites. These five in¬ 
teractions are treated on a mean-field level, whereas the 
antisymmetric DM interactions are neglected for simplicity, 
as they would affect only the lowest frequency range of the 
spectrum in the vicinity of the zone center. This same model 
has been discussed in detail in Ref. 13, for details see also 
Supplementary Information. Having a complete picture of 
magnetic excitations in the entire Brillouin zone up to high 
energies, we can refine the previously established interac¬ 
tion parameters determined from ab-initio calculations^^ 
and ESR experiments^^. 

The values for weak interaction parameters were found 
to be = 27 K, J™ = -50 K and = 45 K, as defined 

W W U..U 


earlier by fitting the magnetization data^l These parame¬ 
ters, in fact, provide low-lying modes in agreement with the 
TOF data as shown in Fig. 5a. Strong coupling constants 
within the tetrahedra, and J™, are mainly responsi¬ 
ble for the positions of high-energy modes and govern the 
intra-tetrahedron excitations. We keep = 145 K as it 
was determined from fitting the ESR spectrum^^. Modifying 
J™ to -170K, we can reproduce the high-energy modes 
in excellent agreement with experiment, as seen in Fig. 5a, 
where the result of the spin-dynamical calculations is shown 
with dotted lines (see also Supplementary Information). 

To fully test our model, we also calculated the scattering 
cross section. At zero temperature this corresponds to 

( 2 ) 

a, (5 

where a,(5 = {x,y,z}, = k(^/\k\, and S“^(k,co) is the 

dynamical spin structure factor. Here we used the origi¬ 
nal crystallographic unit cell and the corresponding cubic 
Brillouin zone, taking the exact positions of each copper 
atom within the unit cell into account. The theoretically 
established scattering intensity plot (artificially broadened 
to model the experimental resolution) is shown in Fig. 5b, 
demonstrating strikingly good agreement with the INS data. 

Conclusions. We have presented the complete overview 
of spin excitations in Cu 20 Se 03 throughout the entire Bril¬ 
louin zone and over a broad energy range. The data reveal 
low-energy ferromagnetic Goldstone modes and a higher- 
energy band of multiple intertwined dispersive spin-wave 
excitations that are separated by an extensive energy gap. 
All the observed features can be excellently described by the 
previously developed theoretical model of interacting CU 4 
tetrahedra, given that one of the strong interaction parame¬ 
ters, J™, is set to a new value of -170 K. The complete list 
of the dominant parameters for the magnetic Hamiltonian 
determined here is now able to describe all the available 
experimental results (magnetization, ESR and INS data) si¬ 
multaneously. Our model can serve as a starting point for 
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Figure 5 | Symmetrized TOF data along high-symmetry directions, a, Energy-momentum profile of the symmetrized low-temperature 
TOF data (T = 4.5 K) along a polygonal path involving most of the high-symmetry directions in Q-space. The data within every segment 
have been averaged with all equivalent cuts in the same Brillouin zone as described in the Methods. The peak positions extracted from the 
fits to the thermal- and cold-neutron TAS data between the r'(222) and X'(220) points (Fig. 2a, b) are shown here with white and green 
data points. The vertical bars marked as a-d on the right-hand side of the panel show energy integration ranges of the corresponding 
panels in Fig. 4. Dotted lines show the calculated dispersions, their transparency reflecting the dynamical structure factors (calculated 
intensities) as described in the text, b. Calculated scattering intensity presented as a color map that was obtained by broadening the 
results of our spin-dynamical calculations with a Gaussian function meant to model the experimental resolution, c. Exchange paths in 
Cu 20 Se 03 . Top: first-neighbor ferromagnetic (blue) and antiferromagnetic (red) couplings. Weak bonds are indicated with dashed lines. 
Bottom: further-neighbor interaction realized on the diagonals of hexagons formed by alternating Cu^ and Cu 2 sites. 


more elaborate low-energy theories that would be able to ex¬ 
plain the complex magnetic phase diagram of Cu 20 Se 03 , in¬ 
cluding the helimagnetic order and skyrmion-lattice phases, 
which still represents a challenge of high current interest in 
solid-state physics. 

Methods. 

Sample preparation. The compound Cu 20 Se 03 has first been synthesized 
by a reaction of CuO (Alfa Aesar 99.995%) and Se 02 (Alfa Aesar 99.999%) 
at 300 °C (2 days) and 600 °C (7 days) in evacuated fused silica tubes. 
Starting from this microcrystalline powder, single crystals of Cu 20 Se 03 
were then grown by chemical transport reaction using NH 4 CI as a trans¬ 
port additive, which decomposes in the vapour phase into ammonia and 
the transport agent HCl. The reaction was performed in a temperature 
gradient from 575 °C (source) to 460 °C (sink), and a transport additive 
concentration of 1 mg/cm^ NH 4 CI (Alfa Aesar 99.999%)^®. The resulting 
single crystals with typical sizes of 30-60mm^ were selectively character¬ 
ized by magnetization, dilatometry, and x-ray diffraction measurements, 
showing good crystallinity and reproducible behavior of physical properties. 


The crystals were then coaligned into a single mosaic sample for neutron¬ 
scattering measurements with a total mass of approximately 2.5 g using a 
digital X-ray Laue camera. 

Triple-axis INS measurements. The measurements presented in Fig. 2 
have been performed at the PUMA thermal-neutron spectrometer with a 
fixed final neutron wave vector kf = 2.662 A“^ and the PANDA cold-neutron 
spectrometer with kf = 1.4A“\ both at the FRM-II research reactor in 
Garching, Germany. A pyrolytic graphite or cold beryllium filter, respec¬ 
tively, was installed after the sample to suppress higher-order scattering 
contamination from the monochromator. The sample was mounted in the 
(HHL) scattering plane and cooled down with a closed-cycle refrigerator 
to the base temperature of ~5 K in both experiments. 

Time-of-flight measurements. TOF data in Figs. 3-5 were collected using 
the Wide Angular Range Chopper Spectrometer (ARCS)^^ at the Spallation 
Neutron Source, Oak Ridge National Laboratory (ORNL), with the incident 
neutron energy set to E[ = 70 meV Here the base temperature of the sample 
was at 4.5 K. The reciprocal space was mapped by performing a full 360° 
rotation of the sample about the vertical ( 110 ) axis in steps of 1 °, counting 
about 20 min per angle. The experimental energy resolution measured as 
the full width at half maximum of the elastic line was ~2.5 meV 
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Data analysis. The TOP data were normalized to a vanadium reference 
measurement for detector efficiency correction, combined, and transformed 
into energy-momentum space using the open-source Matlab -based Horace 
analysis software^^. Due to the high symmetry of the cubic lattice, the 
TOP dataset could be symmetrized in order to improve statistics in the 
data and thereby improve the signal-to-noise ratio considerably. After 
it has been established that the intensity of the INS signal follows the 
symmetry of the large unfolded Brillouin zone introduced in Pig. lb, we 
have assumed this symmetry for the symmetrization. This means that every 
energy-momentum cut in Pig. 3 and every segment of Pig. 5a has been 
averaged with all possible equivalent cuts within the same Brillouin zone 
(cuts from different Brillouin zones with different |Q| were not averaged). 
Por example, the LW or LK segments forming the radius and apothem of 
the hexagonal face of the unfolded Brillouin zone boundary, respectively, 
(Pig. lb) have 48 equivalent instances along 6 nonparallel directions in 
the Brillouin zone that can be averaged. As this kind of data analysis is 
very time consuming and requires significant computational efforts, it still 
awaits to become the standard practice in neutron-scattering research. 
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Supplementary Information to the Letter 
“Magnon spectrum of the helimagnetic insulator Cu20Se02’^ 


Comparison of the theoretical models. 

In Fig. SI, we compare our time-of-flight neutron data with 
theoretical calculations done using the originally published 
exchange parameters from Ref. 13 (top) and the modified 
model with J™ = —170 K (bottom) that demonstrates a sig¬ 
nificantly better agreement with the experimental data, es¬ 
pecially in the regions marked with red arrows. In particular. 


in the high-energy region the old model produced a nearly 
dispersionless intense feature that coincided with a mini¬ 
mum in the observed inelastic scattering intensity, whereas 
in the corrected model this magnon branch is shifted towards 
higher energies. 



r(000) X(020) 1/1/(120) K(| 1 0) r(222) L(111) IA'(120) r(222) X'(220) K(| 1 0) r(000) 

Momentum 



r(000) X(020) 1/1/(120) K(| 10) r(222) L(111) 1/1/(120) r(222) X'(220) K(| 10) r(000) 

Momentum 


Figure SI | The same energy-momentum profile of the symmetrized low-temperature TOF data as shown in Fig. 5a of the paper, compared 
with the theoretically calculated dispersions (dotted lines) for two sets of exchange parameters: a, The original set of parameters from 
Ref. 13; b, the modified set of parameters presented in our current work, where the strong ferromagnetic exchange interaction J™ has 
been modified to -170 K, resulting in a better match with the experimental data. 


In Fig. S2, we additionally compare the experimental line 
profiles at the W points with the corresponding simulated 
profiles obtained from the calculated magnon dispersions 
after resolution broadening (as shown in Fig. 5b in the main 
text). Because the covered energy range for the W(120) 
point is only limited to lower energies, we additionally show 
a cut at the equivalent W^(320) point which provides infor¬ 
mation about the higher-energy part of the spectrum. Both 


points were symmetrized with the same points in equivalent 
Brillouin zones obtained by applying all cyclic permutations 
of {HKL) with possible sign changes to obtain the best signal- 
to-noise ratio. After that both datasets were fitted globally 
with a set of Gaussian peaks to separate the magnetic sig¬ 
nal from the energy-dependent background. This magnetic 
signal is shown as a solid line and compared with the corre¬ 
sponding intensity profile from the theory. 


SI 
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Energy (meV) 

Figure S2 | Comparison of the line profiles at the W(120) and 
1Y(320) points with the corresponding profiles calculated from the 
theoretical model after introducing a finite resolution broadening 
of the magnon peaks. Along with the experimental data (blue 
and red symbols), we additionally show the result of their fitting 
with a set of Gaussian peaks (blue solid line) in order to separate 
them from the energy-dependent background and the model curve 
obtained from the theory (green line). 

Multiboson theory for Cu20Se02 

Having weak and strong interactions suggests that the el¬ 
ementary units are not the bare Cu ions, but the quantum 
mechanically (QM) entangled tetrahedra formed by four 
Cu sites. Therefore, we consider tetrahedra-factorized vari¬ 
ational wave function to describe the ground sate of the 
system and to build upon to introduce the elementary exci¬ 
tations: 

4 

!'?'> = (SI) 

t=l 

where |'0)f is a QM state in the 2^ = 16 dimensional lo¬ 
cal Hilbert space of the strong tetrahedron t. This 16- 
dimensional Hilbert space can be reduced according to the 
C^y local symmetry of the distorted tetrahedron (built of one 
Cui and three Cu 2 sites) as well as the SU(2) rotational sym¬ 
metry. Such classification proves to be useful in labelling our 
states and establishing their mixing without elaborate calcu¬ 
lations. Table SI lists the symmetry classified eigenbasis of 


a strong tetrahedron. 

The ground state is an triplet separated with a large 
energy gap of ~280 K from the rest of the tetrahedron spec¬ 
trum. A finite weak interaction between tetrahedra mixes 
the multiplets, i.e. states with different total spin S. The 
local symmetry, however, remains and the ^ component 
of the spin remains a good quantum number. Thus the A^ 
triplet can only mix with the A^ quintet. Furthermore, if we 
select the = 1 state as the ground state, purely for sake 
of simplicity, the state |1, l,Ai) couples solely to |2, l,Ai) in 
the tetrahedra-factorized picture. The ground state in the 
tetrahedra mean-field case reads 

IV)), = cos ^ |1, l,Ai)t + sin ^|2, l,Ai). (S2) 

The variational parameter a controls the length of local spins 
in the following manner: 

(Sj) = —^ (cosa-h/Ssina) and (S^) = ^ (l — (S^). 

(S3) 

Note that for any a the total magnetization per tetrahedron 
remains (S^ + 3(Sp = 1, i.e. we remain in the one-half 
plateau phase regardless the value of the variational param¬ 
eter a. In the limit a = 0 we recover the decoupled case 
with the local moments (SJ) = 1/4 and (S^) = 5/12. In the 
coupled case for coupling values = 145 K, J™ = —170 K, 
jAF = 27 K, J™ = -50 K, and = 45 K, the variational 
parameter takes the value a = 0.38191 and the local mo¬ 
ments become (Sp = -0.39337 and (S^ = 0.464457. 

To include quantum fluctuations and describe the dynam¬ 
ical properties, we use multiboson theory—a generalizisa- 
tion of the standard spin-wave approach. As a first step, 
we introduce the boson operators ^ that create the 
state (n = 1..16) of a tetrahedron t (t = 1..4). As each 
tetrahedron can only be in one of its 16 states, the bosons 
must satisfy the constraint 

16 

n=l 

The bosonic representation of the spin operator SJ, where 
j denotes a given copper site, is too excessive to show here. 
Nonetheless, some operators acquire a simpler form, such 
as the intra-tetrahedron Hamiltonian and the total magneti¬ 
zation which are diagonal in this basis: 




m=l,0,l 


E=E^,E2 


m=l,0,l 

E=E^,E2 


n,E^l,m,E’ 


(S5) 


Sf (S6) 

The spin raising operator of a tetrahedron t can be written 
as 

^ ^/s{S + 1) — m{m + 1)^5 (S7) 

where the sum is over the 16-dimensional tetrahedron basis, 
and the spin lowering operator is its hermitian conjugate. 


After rewriting our spin operators in terms of boson opera¬ 
tors we transform out the boson that creates the ground 
state to get a Hamiltonian containing only the excitations. 
Then we solve this spin-wave Hamiltonian using the Bogoli- 
ubov transformation to obtain the spectrum. Our ground 
state is expressed in Eqs. SI and S2. Let a[ ^ denote the 
boson that creates IV^), at a tetrahedron t. The orthogonal 
bosons will represent the 15 local excitations on each tetra- 


S2 
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Table SI | The basis of a (strong) tetrahedron classified according to the spin rotational SU(2) and geometrical C^y symmetries. 


irrep 

State 

notation 

al boson 

energy 

Al ® 

5^(1 MiT> +1 UU) +1 fiTif) - 3| nil)) 

^(1 +1 fHiT) +1 ^TTi> -1 ftUT> -1 ftiTf) -1 tTii}) 

^(3i iiTTT) -1 nn) -1 liUT) -1 tm)) 

iiJT 

lhl>Ai 

- 

l,lAi 

<,A. 

p _ 5 tAF _i_ 3 tFM 

E ® 

5^(21 iiitt) -1 nil) -1 UTTf) -1 tUT) -1 nil) + 2 | 1 fTU)) 1 
|(Tnt)-Tm)-iiiUT) + Tm)) j 

\0,0)e 

^0,0,£(1) 
“o,0,E(2, 

T? _ 3 tAF 3 tFM 

^0,£ — 4-'s 4"'s 

Al <g) 

I nil) 

Ki uii) +1 nil) + 1 UTU) +1 iiui)) 

^(l ^iTT) +1 MiT) +1 iiTti) +1 1 iUT) +1 liiTf) +1 nil)) 
Ktttt) + Tm> +1 ftUT) +1 ftTU)) 

1 ftTTT> 

|2,2T 

|2JT 

|2,0)^^ 

|2,1T 

|2,2)^^ 

- 

2,2Ai 

^2,TAi 

2,0Ai 

t 

^2,2Ai 

p _ 3 tAF , 3 tFM 

^2Ai - 4^s + 4^s 


^(TUT) + Tm)-2|iiTU)) 1 

|iJ)e 

- 

iTT(i) 

- 

iTT(2) 


E 0 

^(- 2 | um) +1 iiTiT) +1 UTTf) - 11 iUT) -1 nil) + 2 | 1 fTU)) 1 
KiWT)-Tm> + iftUT)-Tm)) J 

\X0)e 

^1,0,£(1) 
^T0,£(2) 

p _ 1 tAF 3 tFM 

^1,E — 4^s 4^s 


^(-2iifm) + iiinT) + TTn)) \ 

IT1)£ 




hedron. Thus alltogether we have 4 x 15 = 60 excitations. 
To eliminate a[ ^ we use the above constraint and perform 
the following substitution 


a 


t 

IT 



and ^ 



In the linearized multiboson theory we only keep second or¬ 
der boson terms in the resulting Hamiltonian. Therefore, the 
multiboson Hamiltonian, now containing only the excitation 
boson operators, has zero-order boson terms, i.e. constant 
terms linear terms containing single-boson op¬ 
erators, and two-boson operator terms corre¬ 

sponds to the tetrahedral mean-field ground state energy, 
is identically zero when the variational wave function 
minimizes the energy, and the quadratic terms describe 
the boson hopping processes, providing the spectrum. 

In particular, takes the following form: 


^( 2 ) ^ 


1 

2 



where M = M'*’ and N = are 60 x 60 matrices and the 
vector contains the 60 excitations of the unit cell (15 per 
tetrahedron), namely 60 bosons ^ with n = 2,..., 16 and 
t = l,...,4: 




t 


a 


k,16,4 


). 


(SIO) 


The equation of motion reads 

id„^(k) = [a^^(k),^®] = ^(k). (Sll) 

In order to find the linear combinations of bosons that diag¬ 
onalize Jif we need to solve the eigenvalue problem 

The physical eigenvectors (with positive eigenvalues) are 
associated with the states that diagonalize the Hamiltonian 

<k = UT-Al + v^.A_^ (S13) 

Then bosons satisfy the 

commutation relation The Hamiltonian 

becomes diagonal: 

^(2) = 2 + -J • (S14) 


S3 








